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We investigate the possibility of discriminating between Modified Newtonian Dy- 
namics (MOND) and Newtonian gravity with dark matter, by studying the vertical 
dynamics of disk galaxies. We consider models with the same circular velocity in the 
equatorial plane (purely baryonic disks in MOND and the same disks in Newtonian 
gravity embedded in spherical dark matter haloes), and we construct their intrinsic 
and projected kinematical fields by solving the Jeans equations under the assumption 
of a two-integral distribution function. We found that the vertical velocity dispersion 
of deep-MOND disks can be much larger than in the equivalent spherical Newtonian 
models. However, in the more realistic case of high-surface density disks this effect 
is significantly reduced, casting doubts on the possibility of discriminating between 
MOND and Newtonian gravity with dark matter by using current observations. 



ABSTRACT 



Key words: gravitation — 



1 INTRODUCTION 

The flatness of the rotation curves of disk galaxies in their 
external regions has been the primary focus in the last- 

' ing on dark matter (DM) and Modified Newtonian Dy- 
namics (MOND). As a consequence of the formalisation 
of MOND in fundamental physics (Bekenstein 2004; Zlos- 

, nik, Ferreira & Starkman 2007), combined with some dif- 
ficulties of the cold dark matter scenario on galaxy scales 
(e.g. Binney 2004, and references therein), the debate ex- 
tended to other astrophysical contexts. For example, MOND 
models have been recently tested against the Cosmic Back- 
ground Radiation (Skordis et al. 2006, Skordis 2006), grav- 
itational lensing (Chen & Zhao 2006; Zhao et al. 2006; 
Angus et al. 2007; Takahashi & Chiba 2007), star cluster 
and galaxy dynamics (Zhao & Tian 2006; Sanchez- Salcedo, 
Reyes-Iturbide & Hernandez 2006; Haghi, Rahvar & Hasani- 
Zonooz 2006; Nipoti, Londrillo & Ciotti 2007; Scarpa et 
al. 2007; Tiret & Combes 2007), and the solar system 
(Bekenstein & Magueijo 2006, Sanders 2006, Sereno & Jet- 
zer 2006). Due to the surprising ability of MOND to repro- 
duce the kinematics of different systems (e.g., Milgrom 2002; 
Sanders & McGaugh 2002; Bekenstein 2006), it is of obvious 
interest to look for other tests to discriminate between DM 
and MOND. 

In this paper we explore the possibility of differentiat- 
ing DM and MOND by studying the vertical dynamics of 
disk galaxies. While it is often taken for granted that — as 
far as disk kinematics is concerned — the DM and MOND in- 
terpretations are nearly degenerate, in fact this is not true. 
An important issue is the choice of the Newtonian system 
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used for comparison with the MOND results. In the present 
context, the proper comparison is between a baryonic disk in 
MOND, and the same disk, in Newtonian gravity, immersed 
in a DM halo (which for simplicity we assume spherically 
symmetric) , such that the circular velocity in the equatorial 
plane is the same as in the MOND modejj. In particular, we 
compare the kinematical fields of equivalent models under 
the assumption of two-integral distribution function. A work 
in some respect complementary to ours, though with impor- 
tant differences, was done by Milgrom (2001), who studied, 
for a given baryonic disk, the shape of the DM halo of the 
Newtonian model with exactly the same dynamics (not only 
the same rotation curves) as the disk in MOND gravity. Re- 
cently, a test of MOND based on the comparison of disk 
circular velocity and vertical velocity dispersion has been 
proposed also by Stubbs & Garg (20050- 

The paper is organised as follows. In Section [5] we de- 
scribe the general method adopted, and we perform a pre- 
liminary investigation of the problem. In Section[3]we study 
in detail the projected kinematical fields of Miyamoto-Nagai 
and thick exponential disks in deep MOND regime, and we 
compare them with their equivalent Newtonian models with 
DM, for different values of flattening/thickness. We also ad- 



1 Other cases are known in which MOND and Newtonian systems 
are equivalent with respect to some dynamical properties, but 
different with respect to others (e.g. Ciotti & Binney 2004). 

2 The presented test is affected by the unjustified assumption 
that the vertical velocity dispersion depends only on the local 
gravitational field. 



© 2007 RAS 



2 Nipoti et al. 



dress the problem of MOND dynamics in Milky- Way like 
galaxy models. Our results are summarised in Section [5] 



2 THE METHOD 

In the present work we consider MOND in Bekenstein & 
Milgrom's (1984) formulation, in which the Poisson equation 



V 2 <?!> N = 4nGp, 



(1) 

where N is the Newtonian gravitational potential gener- 
ated by the density distribution p, is substituted by the 
non-relativistic field equation 



V ■ 



/' 



Oo 



v<t> 



(2) 



In equation above ||...|| is the standard Euclidean norm, <j> is 
the MOND gravitational potential produced by p, and for 
finite mass systems \7<f> — > for ||x|| — > oo. The interpolating 
function p(t) is not constrained by theory except that it 
must run smoothly from p{i) ~ t at t <C 1 to p(t) ~ 1 at t 
1, with a dividing acceleration scale ao — 1.2 x lCF 10 ms~ 2 . 
In the so-called 'deep MOND regime' (hereafter dMOND), 
describing low-acceleration systems (||V0|| <C ao), p(t) = t 
and equation © reduces to 



V- (||V0||V0) =47rGa 0/ 9. 



(3) 



As well known, equation @ can be combined with equa- 
tion JTJ| so that the MOND g = -V0 and Newtonian 



g 



-V</> N gravitational fields are linked as 



Ms/ao)g = g N + S, 



(4) 



where S is a solenoidal field dependent on the specific p con- 
sidered. It can be proved that the equation above reduces to 
Milgrom's (1983) empirical relation (i.e., S = 0) in case of 
one-dimensional symmetries and in the special case of the 
razor-thin Kuzmin (1956) disk density distribution (Brada 
& Milgrom 1995), but in general one cannot impose S = 0. 
The presented results are derived by solving the field equa- 
tion (O with the numerical potential solver presented in 
Ciotti, Londrillo & Nipoti (2006). 



2.1 Two-integral dynamics 

In the following we compare the dynamics of two disk galaxy 
models with the same circular velocity curve in the equato- 
rial plane 



v!(R) = R- 



dR ' 



(5) 



where 4>t(R,z) is the total potential of the model, and 
(R, z, (fi) are the standard cylindrical coordinates. One model 
is built in MOND for a baryonic disk galaxy described by 
the density distribution p(R,z). The associated Newtonian 
model, which we call equivalent, consists of the same bary- 
onic disk plus a spherical DM halo, whose radial density 
profile is determined by the condition of matching the ro- 
tation curve obtained in the MOND model. In practise, the 
density distribution of the DM halo of the equivalent model 
is given by 



Pdm(t) = 



d 



AnGr 2 dr 



[rv 2 (r) - rv 2 N (r)] , 



(6) 



where v c and v c n{R) are the disk circular velocity in MOND 
and Newtonian gravity, respectively. Note that the positiv- 
ity of pdm is not guaranteed for a generic disk density, so 
we always check the positivity of the equivalent DM halo 
density distribution. 

For a two-integral distribution function, i.e., / = 
f(E,J z ), the Jeans equations for the disk are 



dpa 2 
dz 



= -p- 



djn 

dz 



dpa 2 : p(a 2 -v 2 ) 
dR 



+ 



R 



= ~P- 



d(f)T 
OR 



(7) 



(8) 



where a = a' R = crj, v v — u, v -r u v , u, v 
bar over a symbol indicates its phase-space mean over the 
velocity space (e.g., see Binney & Tremaine 1987, hereafter 
BT). Note that in a MOND model T is the solution of 
equation ([2j for the disk density distribution, while in the 
equivalent Newtonian model 4>t = (f> (R, z) + 0dm (r), with 
4>dm determined from the density profile in equation (]6]). 

The velocity dispersion a 2 (R,z) is obtained from inte- 
gration of equation © with boundary condition pa 2 = for 



pa = 



d<j>T , i 
oz' 



In the equivalent model 



pa = 



oz' 



igDM Z 

dr r 



dz' , 



(9) 



(10) 



where r = \f R 2 + z' 2 and d(f>vM/dr — (v 2 — v 2 N )/r. 

To split v 2 , into streaming motion u v (that for simplicity 
we assume nowhere negative), and azimuthal dispersion a 2 ,, 
we adopt the Satoh (1980) ^-decomposition 



u\ = k 2 (v 2 
and 



a 2 , = a 2 + (1 - fc 2 )(v! 



2\ 



(11) 



(12) 



with < k < 1; this procedure can be applied only when 
v^ — a 2 > everywhere. For k — no ordered motions are 
present, and the velocity dispersion tensor is maximally tan- 
gentially anisotropic, while in the isotropic rotator (k — 1) 
the galaxy flattening is due to azimuthal streaming velocity. 
In principle, by allowing for k = k(R, z), even more rotation- 
ally supported models can be constructed, up to the maxi- 
mum rotation case considered in Ciotti & Pellegrini (1996), 
where a'i = everywhere. 

The explicit projection formulae for the kinematical 
fields of axisymmetric two-integral systems viewed along a 
generic direction n of the line of sight (los) can be found 
elsewhere (e.g., see Lanzoni & Ciotti 2003, Ciotti & Bertin 
2005, Riciputi et al. 2005): here we just recall that at any 
given place in the galaxy the los component of the stream- 
ing velocity field is tt n = (v, n) = UiUi, where (, ) is the 
standard inner product, the repeated index convention is 
applied, and Cartesian coordinates are assumed. The analo- 
gous quantity associated with the velocity dispersion tensor 
is (j 2 = (v — u) 2 = aijUiTij, and the corresponding (mass- 
weighted) kinematical projected fields are obtained by inte- 
gration along the los as 



pu n dl, 



(13) 
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Figure 1. Density profile of spherical DM haloes of Newto- 
nian models equivalent to dMOND MN disks (left) and thick 
exponential disks (right). The dotted lines in the left panel are 
given by the distribution i ]24[ ). obtained in the S = approxi- 
mation. Density is in units of (M/4nb 3 ) ^Jaob 2 /GM (left) and 
(M/4mRlz )y/aoRyGM (right). 



Figure 2. Ratio of the vertical accelerations of dMOND (g z disk) 
and equivalent Newtonian (g^ DM ) models as a function of the 
height z above the disk plane, at fixed cylindrical radius. 



where simp — y/R = y/\/y 2 + x 2 . Note that, after projec- 
tion, the y coordinate can been identified with R. 



and 



E<7„ 



pu a dl, 



paldl, 



(14) 



(15) 



where E = J_ p dl is the disk density projected along n. 
In general <r p is not the velocity dispersion measured in ob- 
servations: in the presence of a non-zero projected velocity 
field v p , this quantity is 



EaL = P « v > n > - v p ) 2 dl = E x (4 + V 2 - v 2 p ) . (16) 

J — oc 

For simplicity, in this paper we focus on the limit cases 
of disk observed face-on and edge-on. The face-on projected 
velocity dispersion is given by 

poo f OC C\l 

Ecr 2 = 2 / pa 2 dz = 2 p-^zdz, (17) 
Jo Jo " z 

where E(i?) = 2 J"°° p{R, z)dz is the face-on disk density. In 
the edge-on case, assuming the los directed along the x axis, 
the los direction in the natural coordinate system is given 
D J0 n — (1) 0, 0) and the projection plane is (y, z), so that 



u n = — u v sm ip, 
and 



o"n = + (1 — fc 2 )(v^> — <r' a ) sin^ <p 



(18) 
(19) 



3 The los vector points towards the observer, and so positive ve- 
locities correspond to a blue-shift. 



2.2 A preliminary analysis 

As an introductory exercise we compare, as a function of 
disk flattening, the vertical force field near the equatorial 
plane of a Miyamoto & Nagai (1975, hereafter MN) disk 
in dMOND regime and its equivalent Newtonian model. In 
order to carry out the calculations analytically we assume 
S = 0, while in Section [3] the exact MOND gravitational 
field with the appropriate S field is used. The MN density 
distribution with scale length b and flattening parameter 
s = a/b is 



p(R,z) = 



M sR 2 + (s + 3C)Q + C) 2 
4tt6 3 [£a + ( s + £)2]5/2£3 



(20) 



where M is the total disk mass, Q = \/l + z 2 , z = z/b, and 
R = R/b. The corresponding Newtonian potential is (BT) 

GM 1 



<i>(R,z) 



b ^R 2 + (s + C) 2 ' 



(21) 



For a = the MN density distribution is a Plummer (1911) 
sphere with scale radius b, while for 6 = one obtains the 
Kuzmin (1956) disk. The Newtonian disk circular velocity 
is given by equations <(Sj and (|21|l as 



Vcn(R) 



GM 



R 2 



[R 2 + {s + l) 2 ] 3 / 2 



(22) 



In the S = approximation the dMOND circular velocity is 
obtained from equation @ as 



2 . m VGMaTiR^ 2 
v c (R) — — — , 

[R 2 + (S + l)2]3/4' 

and according to equation 



(23) 



©, pom oc A(r) 



y GM/aob 2 B(f), where A and B are two dimensionless 
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Figure 3. Radial trend of velocity ratios between dMOND and equivalent Newtonian models (quantities with label N) in the fully 
isotropic case. From top to bottom: 2 = vertical velocity dispersion l[9j, face-on projected velocity dispersion dlTI . z = edge-on 
projected streaming velocity ( 1131 . and z = edge-on projected velocity dispersion 116> . 



functions, and r = r/b. In the dMOND limit yjGM/a b 2 = 
and so 



Pdm(V) 



M 



an 



b 2 2f 2 + 5(s + iy 



4tt6 3 V GM 2\f?\r 2 + (s + l) 2 ] 7 / 4 ' 



(24) 



which is positive everywhere and has the expected r~ 2 be- 
haviour for r — *■ co. Note that equation (|24l) implies that in 
the equivalent Newtonian model the MN disk is just a tracer 
of the DM gravitational field. 

The obtained formulae are useful because they allow to 
evaluate the vertical gravitational field near the disk plane 
in the two models. In the S = approximation the dMOND 
vertical field is 



Sz.disk 



ao 

W\ 



■9z , 



(25) 



where g^ is the vertical Newtonian gravitational field of the 
MN disk, and the vertical gravitational field of the equivalent 
Newtonian DM halo is 



N 

<?z,DM — 



dcj>DM(r) _ zv 2 (r) 



dz r 2 
Simple algebra then shows that for R 



(26) 



oo and z — > 0, 



Sz,DM 



1 



Sz,disk 1+s 



+ 



s(5 + s) 



2(1 



(27) 



' 4(l + s)R 2 _ 

i.e., the vertical field is stronger in dMOND than in DM 
models, and the discrepancy between the two cases is larger 
for flatter models (i.e. for larger values of s). Thus, we ex- 
pect a 2 in the disk to be different in the two models. In 
turn, equations JH} and (fTTjl imply that also the streaming 



velocity field in the disk equatorial plane is different in the 
two models, even though the circular velocity is the same by 
construction. 



3 RESULTS 

We now obtain the kinematical fields (for simplicity in the 
fully isotropic case) of a few disk galaxy models by comput- 
ing the exact MOND and Newtonian potentials with the nu- 
merical code presented in Ciotti et al. (2006) . The potential 
solver is based on a spherical grid of coordinates (r, •&, tp), 
with N r x N# x iV v points, uniformly spaced in (arctanr, #, 
ip). To compute the intrinsic and projected velocity disper- 
sions, we interpolate the potential from the spherical grid to 
a cylindrical grid with Nr x N z points, uniform in (arctan R, 
arctan z), on which we evaluate the integrals (|9} and (|13H15[) . 
The numerical integration routines are verified by deriving 
the intrinsic and projected velocity dispersions of MN disks 
in Newtonian gravity, and comparing the results with the an- 
alytical expressions reported in Ciotti & Pellegrini (1996). 
For instance, with iV r = N$ = Nr = N z = 128 the analyti- 
cal results are reproduced with relative errors ~ 0.1 %. 



3.1 Miyamoto-Nagai disks 

We consider two MN galaxy models with flattening param- 
eter s = 1 (i.e., an almost spherical system), and s = 5. 
Their exact dMOND fields are computed by solving numer- 
ically equation @. The density distributions of the spheri- 
cal DM haloes of the equivalent Newtonian models are then 
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Figure 4. Density distribution in the equatorial plane for the 
baryonic and dark components of the considered Milky-Way like 
galaxy models (top), and the corresponding rotation curves (bot- 
tom). 



obtained from equation (JS} with u c n = 0, as appropriate 
for the dMOND case, in which the disk is just a tracer of 
the halo potential. As expected, the resulting halo profiles 
(Fig. \T\ left panel; solid and dashed lines) decrease as r~ 2 
in the outer regions, matching the analytical formula (|24[) . 
obtained in the S = approximation (dotted lines). How- 
ever, for r<b the exact halo density profiles are consider- 
ably natter than those predicted by equation ([24|l . having 
inner logarithmic slope d log p^m/d log r> — 1/2. The differ- 
ence (which is larger in the flatter disk) indicates that the 
contribution of the field S is non-negligible in the central 
regions, as also found by Ciotti et al. (2006) for triaxial and 
axisymmetric Hernquist (1990) models in dMOND regime. 
However, in the outer regions the agreement is perfect. 

In Fig. [2] (left panel), we plot the ratio <?^DM/5z,disk as 
a function of z at the representative radius R — b. Again, 
the results nicely confirm the asymptotic estimate of equa- 
tion (|27[) when z = 0: the effect is bigger in the flatter disk, 
and just above the disk plane <?z,disk is significantly stronger 
than 5™ DM . We also note that for z 2> b the vertical New- 
tonian field exceeds the dMOND one, though slightly. As a 
consequence, also the intrinsic and projected velocity disper- 
sions are different in the equivalent Newtonian and dMOND 
MN disks, while the circular velocities in the equatorial plane 
are the same by construction. The midplane vertical veloc- 
ity dispersion u(R, 0) and the projected face-on velocity dis- 
persion crios are larger in the dMOND systems than in the 
equivalent Newtonian systems at all radii R < 5b, with the 
larger discrepancies (up to a factor of 2) in the natter disk 
(Fig. [3j left column). The situation is more complicate in 
the edge-on projection: the projected streaming velocity v p 
in the Newtonian models are larger than in the dMOND 
cases, while the edge-on projected velocity dispersion <ti os in 
the galactic plane is higher in dMOND than in Newtonian 
models. 



Figure 5. Vertical force as a function of z, at R = 8kpc, for the 
considered Milky- Way like galaxy models. 

3.2 Thick exponential disks 

We now investigate the more realistic thick exponential disk 

^^-4^ exp (-|) sech2 (i)' (28) 

where M is the total disk mass, and Ro and zo are character- 
istic scale- lengths. As for the MN case, we consider here two 
dMOND models (and so in the equivalent Newtonian model 
the disk does not contribute to the potential) : a thicker disk 
with Ro — 3zo, and a thinner disk with Ro = 15zo. 

As in MN disks, the spherical DM haloes density pro- 
files of the equivalent Newtonian models are rather flat in the 
central regions and oc r~ 2 at large radii (Fig.[T] right panel). 
Interestingly, these DM haloes "predicted" by MOND for 
low-surface brightness disks are not characterised by the 
steep central cusps expected in the context of cold dark mat- 
ter. The dMOND vertical force near the disk is significantly 
stronger than in the equivalent Newtonian models (Fig. [3] 
right panel): at R — Ro in the thinner model the dMOND 
vertical force ratio is ~ 10 just above the plane. This differ- 
ence reflects in the velocity fields (Fig. [31 right panels), which 
behave qualitatively as the MN models, but the discrepan- 
cies are larger: for example, the ratio of intrinsic velocity 
dispersion is as high as ~ 4 in the thinner model. 

3.3 Milky- Way like galaxies 

The previous analysis restricted to models in the dMOND 
regime, in order to estimate the largest possible differences 
in the kinematical fields of equivalent disk galaxies. We now 
focus on a Milky- Way like galaxy model, with the purpose of 
quantifying the MOND kinematical effects in realistic high 
surface-brightness galaxies. We will not enter the challeng- 
ing problem of finding the best-fitting MOND model of the 
Milky Way (which is beyond the aim of the present paper; 
see Famaey & Binney 2005); instead, we apply our method 
to one of the currently available mass models of the Milky 
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Figure 6. Left (from top to bottom): z = vertical velocity dispersion, face-on projected velocity dispersion, z = edge-on projected 
streaming velocity and 2 = edge-on projected velocity dispersion for model DBMOND (solid lines) and its equivalent Newtonian model 
DBEQ (dashed lines), in units of kms -1 . Right: percentage absolute differences between the corresponding MOND and Newtonian 
velocity fields in the left panels. Thin and thick lines refer to thin-disk and thick-disk stars, respectively. 



Way that have been derived in the context of DM. In par- 
ticular, we consider Model 1 of Dehnen & Binney (1998; 
hereafter model DB), which consists of a truncated oblate 
spheroidal power-law bulge, a stellar disk made of a thin and 
a thick component, a gaseous exponential disk with a cen- 
tral hole, and finally an oblate spheroidal DM distribution. 
Model DB represents a good Newtonian dynamical model of 
the solar neighbourhood and the Galaxy, and can be classi- 
fied as a disk dominated model, because in Newtonian grav- 
ity the disk provides the main contribution to the rotation 
curve in the central regions. 

As we solve now equation ([2} , we need to specify the in- 
terpolating function /i. Following Famaey & Binney (2005), 
we fix ao ~ 1.2 x 10™ 10 ms~ 2 and 

/*(*) = YTt- (29) 

This function gives a better fit to the terminal velocity 
curve of the Milky Way than the commonly adopted func- 
tion /j,(t) — t/y/1 + t 2 and fits extremely well the rotation 
curve of NGC3198, also being consistent with Bekenstein's 
(2004) TeVeS theory (Zhao & Famaey 2006; see also Famaey 
et al. 2007). We consider two models with the same bary- 
onic distribution as model DB: model DBMOND (a purely 
baryonic MOND model), and model DBEQ (the equiva- 
lent Newtonian model with spherical DM halo). Clearly, 
the DM halo density profile in model DBEQ is not iden- 
tical to that of model DB: Fig. [4] (top panel) shows that in 
the equatorial plane the two distributions are very similar 
in the radial range R ~ 5 — 15 kpc, while the halo profile 



of DBEQ is steeper than that of DB at small radii. As a 
consequence, the rotation curve of models DBMOND and 
DBEQ (which is the same by construction) differs from the 
rotation curve of DB, being systematically higher at small 
radii and lower at large radii (Fig. 3] bottom). However, 
the differences in the radial range — 15 kpc are within 
~ 5%, and u c (8kpc) = 226 kms" 1 in model DBMOND (and 
DBEQ), consistent with the observational constraint at the 
solar neighbourhood 220 ± 15 kms -1 (solid symbol and ver- 
tical bar in the diagram; e.g., BT). 

In contrast with the dMOND cases discussed in Sec- 
tions [37T] and [331 now the contribution of the baryonic com- 
ponent to the gravitational field is important also in the 
Newtonian model, and this results in smaller differences be- 
tween the MOND and Newtonian equivalent models. For 
instance, considering the vertical force as a function of z 
at R = 8 kpc (Fig. |3J, we find that the MOND force is 
stronger than the Newtonian one, with a difference of ~ 30% 
at z ~1 kpc. The solid dot marks an observational estimate 
(and associated uncertainty) of the vertical force 1.1 kpc 
above the plane at the solar radius (Kuijken & Gilmore 1989, 
1991). Remarkably, the difference between the vertical forces 
at z = 1.1 kpc predicted by the MOND and Newtonian mod- 
els is larger than the observational error. So, taken at face 
value, this result would favour DM rather than MOND mod- 
els, though we stress again that we did not attempt to build 
the best-fit MOND model of the Milky Way, and we can- 
not exclude the possibility of finding a MOND model of 
the Milky Way satisfying this and other observational con- 
straints (see also Famaey & Binney 2005). In any case, this 
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Figure 7. Isodensity contours in the meridional plane of the fully equivalent DM haloes associated with dMOND MN (top panels) and 
thick exponential disks (bottom panels). Density increases towards the equatorial plane. 



result shows that in principle it could be possible to try to 
falsify MOND by using the observational constraints on the 
vertical force above the plane. 

Following the treatment of Sections 13.11 and 13.21 we 
computed the kinematical fields of models DBMOND and 
DBEQ, again under the same assumptions of two-integral 
distribution functiorQ and full isotropy. For both models we 
computed the intrinsic and projected velocity fields of the 
thin-disk (thin curves) and thick-disk (thick curves) stars 
(Fig. [6] left panels), finding that the discrepancy between 
the two models is up to 20-25% for intrinsic and face-on 
velocity dispersion, while just few per cent in case of edge- 
on projection (Fig. [6] right panels). These results show that 
the differences in the vertical kinematical fields in a realis- 
tic high-surface density galaxy are quite small, and it is not 
obvious that one can discriminate between MOND and New- 
tonian gravity (plus DM) by using the information currently 
derived by the observations. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we quantified the differences in the vertical 
force and kinematical fields between MOND disk galaxy 
models and equivalent (i.e., having identical circular veloc- 
ity in the equatorial plane) Newtonian models with spherical 



4 Because of this assumption our results on the velocity fields 
cannot be intended, strictly speaking, to represent the specific 
case of the Milky Way, whose distribution function is not two- 
integral, being > a\ in the solar neighbourhood. 



DM haloes. We showed that, in principle, MOND and New- 
tonian gravity with DM can be differentiated with accurate 
measurements of the vertical force near the disk and of the 
projected kinematical fields. In particular, MOND models 
have stronger vertical force near the plane and higher ver- 
tical velocity dispersion than equivalent Newtonian models 
with DM. Our results are not inconsistent with the finding of 
Read & Moore (2005) that satellite orbits in DM and MOND 
behave very similarly, because far from the disk plane the 
force fields predicted by the two theories are very similar. 

From the observational point of view, the robust dis- 
crimination between MOND and Newtonian gravity with 
DM is still challenging due to several factors: (i) the 
strongest effects are expected in low-acceleration (low- 
surface-brightness) systems, in which measurements of the 
velocity fields are difficult; (ii) there are often several stel- 
lar populations with different scale-lengths and heights, and 
the mass-to-light of the stellar component is tunable; (iii) 
most disks have gaseous components whose density has non- 
smooth features which could make differences in any local 
volume of kpc scale; (iv) the models should be compared to 
data with noise so our comparisons of MOND and DM mod- 
els of exactly same circular velocity curves remains a theo- 
retical exercise; (v) DM haloes are not necessarily spherically 
symmetric. 

With regard to point (v), we briefly discuss how our re- 
sults would change by relaxing the assumption of an equiva- 
lent spherical DM halo. Following Milgrom (2001), we start 
by constructing, for a given disk distribution p, the "fully 
equivalent" DM density distribution as pdm = V 2 <^/47rG— p, 
where cfi is the MOND potential; as well known, pdm can be 
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negative in some region of space (Milgrom 1986; Ciotti et 
al. 2006). By construction, the fully equivalent Newtonian 
model has gravitational field identical to the MOND disk 
field in all the space, not just the same midplane circular 
velocity (as the equivalent model with spherical halo). In- 
terestingly, we found that the density distributions poM of 
the fully equivalent DM haloes of the dMOND models pre- 
sented in Sections l3.1l and l3.2l are nowhere negative (in these 
cases pDM = V 2 </>/47rG, because the baryonic disk is just a 
tracer in dMOND-equivalent models). Figure [7] shows that 
Pdm presents a significant disk structure (more pronounced 
in the case of more flattened baryonic disks). This result 
clearly indicates that it is not possible to obtain a Newto- 
nian model with spheroidal DM distribution strictly equiva- 
lent to a MOND disk model (see also Milgrom 2001). How- 
ever, we found that the discrepancy in the vertical dynamics 
in the galactic plane between MOND and DM models can 
be significantly reduced considering spheroidal DM haloes: 
for instance, a dMOND MN disk with s = 5 and an oblate 
spheroidal halo with axis ratio ~ 0.25 — 0.3 (and the same 
asymptotic circular velocity as the baryonic dMOND disk) 
would produce very similar vertical velocity dispersion and 
rotational velocity in the galactic plane (excluding the very 
central regions). 

For all these reasons, it appears difficult to differenti- 
ate between MOND and DM using the currently available 
observational measures of stellar kinematics of external disk 
galaxies (e.g., Bottema 1993; Kregel & van der Kruit 2005, 
and references therein). The most promising application of 
our method seems to be the case of the Milky Way: in partic- 
ular, for the specific baryonic mass distribution considered in 
this paper the difference between the MOND and DM ver- 
tical forces above the plane is larger than the uncertainty 
on the observational estimates. The estimates of both the 
baryonic mass distribution and the vertical force above the 
plane are expected to be greatly improved in the near future 
thanks to data from the GAIA mission. 
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